"""
Comprehensive test suite for advanced PDE implementations

This test suite validates the advanced PDE types and compares results
with py-pde implementations where possible.
"""

using Test
using PDEJulia
using PDEJulia.PDEs
using PDEJulia.Fields
using PDEJulia.Grids

# Test configuration
const TEST_TOLERANCE = 1e-10
const INTEGRATION_TIME = 1.0
const GRID_SIZE = 32

@testset "Advanced PDE Implementation Tests" begin
    
    @testset "Cahn-Hilliard PDE Tests" begin
        @testset "Basic Functionality" begin
            # Create test grid and initial condition
            grid = UnitGrid([GRID_SIZE])
            x = range(0, 1, length=GRID_SIZE)
            u0_data = 0.1 * cos.(2π * x) + 0.05 * sin.(4π * x)
            u0 = ScalarField(u0_data, grid)
            
            # Create Cahn-Hilliard PDE
            pde = CahnHilliardPDE(interface_width=0.1, mobility=1.0)
            
            @test isa(pde, CahnHilliardPDE)
            @test pde.interface_width == 0.1
            @test pde.mobility == 1.0
            
            # Test evolution rate computation
            dudt = evolution_rate(pde, u0, 0.0)
            @test isa(dudt, ScalarField)
            @test size(data(dudt)) == size(data(u0))
            
            # Test mass conservation (∫dudt dV = 0 for Cahn-Hilliard)
            mass_change = sum(data(dudt))
            @test abs(mass_change) < TEST_TOLERANCE
        end
        
        @testset "Energy Functional" begin
            grid = UnitGrid([GRID_SIZE])
            u0 = ScalarField(0.5 * ones(GRID_SIZE), grid)
            pde = CahnHilliardPDE(interface_width=1.0)
            
            energy = free_energy(pde, u0)
            @test energy >= 0.0
            
            # Energy should decrease with time evolution
            dudt = evolution_rate(pde, u0, 0.0)
            u1 = ScalarField(data(u0) + 0.01 * data(dudt), grid)
            energy1 = free_energy(pde, u1)
            @test energy1 <= energy  # Energy non-increasing
        end
        
        @testset "Conservation Properties" begin
            grid = UnitGrid([GRID_SIZE])
            u0_data = randn(GRID_SIZE) * 0.1
            u0 = ScalarField(u0_data, grid)
            pde = CahnHilliardPDE()
            
            @test is_mass_conserved(pde) == true
            
            initial_mass = sum(data(u0))
            dudt = evolution_rate(pde, u0, 0.0)
            
            # Mass conservation: ∫(∂u/∂t) dV = 0
            mass_rate = sum(data(dudt))
            @test abs(mass_rate) < TEST_TOLERANCE
        end
    end
    
    @testset "KPZ Interface PDE Tests" begin
        @testset "Basic Functionality" begin
            grid = UnitGrid([GRID_SIZE])
            h0_data = 0.1 * randn(GRID_SIZE)
            h0 = ScalarField(h0_data, grid)
            
            pde = KPZInterfacePDE(nu=0.5, lambda=1.0, noise=0.0)
            
            @test isa(pde, KPZInterfacePDE)
            @test pde.nu == 0.5
            @test pde.lambda == 1.0
            @test has_noise(pde) == false
            
            dhdt = evolution_rate(pde, h0, 0.0)
            @test isa(dhdt, ScalarField)
        end
        
        @testset "Nonlinear Terms" begin
            grid = UnitGrid([GRID_SIZE])
            x = range(0, 1, length=GRID_SIZE)
            h0_data = sin.(2π * x)
            h0 = ScalarField(h0_data, grid)
            
            pde = KPZInterfacePDE(nu=0.0, lambda=1.0, noise=0.0)  # Pure nonlinear
            dhdt = evolution_rate(pde, h0, 0.0)
            
            # Should have nonlinear gradient contribution
            @test !all(data(dhdt) .== 0)
            @test is_nonlinear(pde) == true
        end
        
        @testset "Stochastic Terms" begin
            grid = UnitGrid([GRID_SIZE])
            h0 = ScalarField(zeros(GRID_SIZE), grid)
            
            pde_noise = KPZInterfacePDE(nu=0.5, lambda=0.0, noise=1.0)
            @test has_noise(pde_noise) == true
            
            # Stochastic evolution should be different each time
            dhdt1 = evolution_rate(pde_noise, h0, 0.0)
            dhdt2 = evolution_rate(pde_noise, h0, 0.0)
            
            # Should be different due to noise (probabilistically)
            difference = norm(data(dhdt1) - data(dhdt2))
            @test difference > 0.1  # Should be significantly different
        end
    end
    
    @testset "Kuramoto-Sivashinsky PDE Tests" begin
        @testset "Basic Functionality" begin
            grid = UnitGrid([GRID_SIZE])
            u0_data = 0.1 * cos.(2π * range(0, 1, length=GRID_SIZE))
            u0 = ScalarField(u0_data, grid)
            
            pde = KuramotoSivashinskyPDE(nu=1.0)
            
            @test isa(pde, KuramotoSivashinskyPDE)
            @test is_chaotic(pde) == true
            @test is_stiff(pde) == true
            
            dudt = evolution_rate(pde, u0, 0.0)
            @test isa(dudt, ScalarField)
        end
        
        @testset "Energy Dissipation" begin
            grid = UnitGrid([GRID_SIZE])
            u0_data = randn(GRID_SIZE) * 0.1
            u0 = ScalarField(u0_data, grid)
            pde = KuramotoSivashinskyPDE()
            
            energy0 = get_energy(pde, u0)
            @test energy0 >= 0.0
            
            # Test dissipation rate computation
            dissipation = get_dissipation_rate(pde, u0)
            @test dissipation >= 0.0  # Should dissipate energy
        end
        
        @testset "Time Scale Estimation" begin
            grid = UnitGrid([GRID_SIZE])
            pde = KuramotoSivashinskyPDE(nu=1.0)
            
            grid_spacing = 1.0 / GRID_SIZE
            time_scale = get_time_scale(pde, grid_spacing)
            @test time_scale > 0.0
            @test time_scale ≈ grid_spacing^4  # Fourth-order scaling
        end
    end
    
    @testset "Swift-Hohenberg PDE Tests" begin
        @testset "Pattern Formation" begin
            grid = UnitGrid([GRID_SIZE])
            u0_data = 0.01 * randn(GRID_SIZE)
            u0 = ScalarField(u0_data, grid)
            
            # Supercritical case (ε > 0)
            pde_super = SwiftHohenbergPDE(rate=0.1, kc2=1.0)
            @test is_pattern_forming(pde_super) == true
            
            # Subcritical case (ε < 0)
            pde_sub = SwiftHohenbergPDE(rate=-0.1, kc2=1.0)
            @test is_pattern_forming(pde_sub) == false
        end
        
        @testset "Critical Wavelength" begin
            pde = SwiftHohenbergPDE(kc2=4.0)  # kc = 2
            wavelength = get_critical_wavelength(pde)
            @test wavelength ≈ π  # λc = 2π/kc = 2π/2 = π
        end
        
        @testset "Linear Growth Rate" begin
            pde = SwiftHohenbergPDE(rate=0.1, kc2=1.0)
            
            # At critical wave number, should have maximum growth
            k_critical = 1.0  # kc
            growth_rate = get_linear_growth_rate(pde, k_critical)
            @test growth_rate ≈ 0.1  # Should equal rate parameter
            
            # At k = 0, should have negative growth
            growth_zero = get_linear_growth_rate(pde, 0.0)
            @test growth_zero < 0
        end
        
        @testset "Energy Functional" begin
            grid = UnitGrid([GRID_SIZE])
            u0 = ScalarField(0.1 * ones(GRID_SIZE), grid)
            pde = SwiftHohenbergPDE(rate=0.1)
            
            energy = free_energy(pde, u0)
            @test isfinite(energy)
        end
    end
    
    @testset "Wave PDE Tests" begin
        @testset "Energy Conservation" begin
            grid = UnitGrid([GRID_SIZE])
            x = range(0, 1, length=GRID_SIZE)
            u0_data = sin.(2π * x)
            v0_data = zeros(GRID_SIZE)
            
            u0 = ScalarField(u0_data, grid)
            v0 = ScalarField(v0_data, grid)
            state0 = FieldCollection([u0, v0])
            
            pde = WavePDE(speed=1.0)
            
            @test is_hyperbolic(pde) == true
            @test requires_cfl_stability(pde) == true
            
            initial_energy = wave_energy(pde, state0)
            @test initial_energy > 0.0
            
            # Evolution should conserve energy (for linear wave equation)
            dstatedt = evolution_rate(pde, state0, 0.0)
            @test length(dstatedt) == 2  # Should return [dudt, dvdt]
        end
        
        @testset "CFL Condition" begin
            pde = WavePDE(speed=2.0)
            grid_spacing = 0.1
            dt = 0.04
            
            cfl, is_stable = cfl_condition(pde, grid_spacing, dt)
            @test cfl ≈ 0.8  # 2.0 * 0.04 / 0.1 = 0.8
            @test is_stable == true  # CFL < 1
        end
    end
    
    @testset "Multi-field PDE Tests" begin
        @testset "Reaction-Diffusion System" begin
            # Test FitzHugh-Nagumo system
            reaction_func = FitzHughNagumoReaction()
            pde = ReactionDiffusionPDE(0.1, 0.01, reaction_func)
            
            grid = UnitGrid([GRID_SIZE])
            u0 = ScalarField(0.1 * ones(GRID_SIZE), grid)
            v0 = ScalarField(0.05 * ones(GRID_SIZE), grid)
            state0 = FieldCollection([u0, v0])
            
            dstatedt = evolution_rate(pde, state0, 0.0)
            @test length(dstatedt) == 2
            @test isa(dstatedt[1], ScalarField)
            @test isa(dstatedt[2], ScalarField)
        end
        
        @testset "Conservation Properties" begin
            # Test Gray-Scott system (mass conserving reactions)
            reaction_func = GrayScottReaction()
            pde = ReactionDiffusionPDE(0.2, 0.1, reaction_func)
            
            grid = UnitGrid([GRID_SIZE])
            u0 = ScalarField(ones(GRID_SIZE), grid)
            v0 = ScalarField(zeros(GRID_SIZE), grid)
            state0 = FieldCollection([u0, v0])
            
            # Total concentration should be conserved in Gray-Scott
            initial_total = get_total_concentration(state0)
            @test initial_total > 0.0
        end
    end
    
    @testset "Stochastic PDE Tests" begin
        @testset "Noise Generation" begin
            grid = UnitGrid([GRID_SIZE])
            u0 = ScalarField(zeros(GRID_SIZE), grid)
            
            # White noise
            white_noise = WhiteNoise(1.0)
            noise_field = generate_noise(white_noise, u0, 0.01)
            @test isa(noise_field, ScalarField)
            @test var(data(noise_field)) > 0.0
        end
        
        @testset "Stochastic Allen-Cahn" begin
            pde = StochasticAllenCahnPDE(temperature=0.1)
            @test isa(pde, StochasticPDE)
            @test has_noise(pde) == true
            
            grid = UnitGrid([GRID_SIZE])
            u0 = ScalarField(zeros(GRID_SIZE), grid)
            
            dudt = evolution_rate(pde, u0, 0.0; dt=0.01)
            @test isa(dudt, ScalarField)
        end
    end
    
    @testset "Expression and Utility Tests" begin
        @testset "PDE Expressions" begin
            pde_allen = AllenCahnPDE(interface_width=2.0, mobility=0.5)
            expr = expression(pde_allen)
            @test isa(expr, String)
            @test contains(expr, "∇²c")
            
            pde_wave = WavePDE(speed=3.0)
            expr_wave = expression(pde_wave)
            @test contains(expr_wave, "v")
            @test contains(expr_wave, "9∇²u")  # speed² = 9
        end
        
        @testset "PDE Classification" begin
            @test is_stiff(CahnHilliardPDE()) == true
            @test is_stiff(DiffusionPDE()) == false
            
            @test is_hyperbolic(WavePDE()) == true
            @test is_hyperbolic(DiffusionPDE()) == false
            
            @test requires_high_accuracy(KuramotoSivashinskyPDE()) == true
        end
    end
    
    @testset "Solver Integration Tests" begin
        @testset "Auto-solver Selection" begin
            # Test that auto-solver works for different PDE types
            grid = UnitGrid([32])
            u0 = ScalarField(0.1 * randn(32), grid)
            
            # Stiff PDE should get implicit solver
            pde_stiff = CahnHilliardPDE()
            solver_config = PDEs.select_optimal_solver(pde_stiff, u0, (0.0, 1.0))
            @test haskey(solver_config, :solver_type)
            
            # Hyperbolic PDE should get explicit solver with CFL
            u_wave = ScalarField(sin.(2π * range(0, 1, length=32)), grid)
            v_wave = ScalarField(zeros(32), grid)
            state_wave = FieldCollection([u_wave, v_wave])
            pde_wave = WavePDE()
            solver_config_wave = PDEs.select_optimal_solver(pde_wave, state_wave, (0.0, 1.0))
            @test haskey(solver_config_wave, :cfl) || haskey(solver_config_wave, :dt)
        end
    end
end

@testset "Comparison with py-pde Reference" begin
    # These tests compare against known analytical solutions or
    # py-pde reference implementations where available
    
    @testset "Diffusion Equation Comparison" begin
        # Test against analytical solution: u(x,t) = exp(-π²Dt) sin(πx)
        grid = UnitGrid([64])
        x = range(0, 1, length=64)
        D = 0.1
        
        # Initial condition: sin(πx)
        u0_data = sin.(π * x)
        u0 = ScalarField(u0_data, grid)
        
        pde = DiffusionPDE(diffusivity=D)
        
        # Analytical solution at small time
        t = 0.01
        analytical = exp(-π^2 * D * t) * sin.(π * x)
        
        # Numerical evolution (simplified - would need actual time stepping)
        dudt = evolution_rate(pde, u0, 0.0)
        u_approx = data(u0) + t * data(dudt)
        
        # Check that we're evolving in the right direction
        @test norm(u_approx - analytical) < norm(data(u0) - analytical)
    end
    
    @testset "Wave Equation Comparison" begin
        # Standing wave solution: u(x,t) = A cos(ωt) sin(kx)
        grid = UnitGrid([64])
        x = range(0, 1, length=64)
        
        # Initial conditions for standing wave
        A, k, c = 1.0, π, 1.0
        ω = c * k
        
        u0_data = A * sin.(k * x)
        v0_data = zeros(64)  # Zero initial velocity
        
        u0 = ScalarField(u0_data, grid)
        v0 = ScalarField(v0_data, grid)
        state0 = FieldCollection([u0, v0])
        
        pde = WavePDE(speed=c)
        
        # Test that energy is initially correct
        initial_energy = wave_energy(pde, state0)
        expected_energy = c^2 * k^2 * A^2 * length(x) / 4  # Simplified
        
        # Should have right order of magnitude
        @test initial_energy > 0.0
        
        # Evolution should maintain the wave character
        dstatedt = evolution_rate(pde, state0, 0.0)
        @test data(dstatedt[1]) ≈ data(v0)  # du/dt = v
    end
end

println("All advanced PDE tests completed!")